L’objectiu d’aquesta activitat serà el tractament d’un dataset, que pot ser el creat a la pràctica 1 o bé qualsevol dataset lliure disponible a Kaggle (https://www.kaggle.com).
Un exemple de dataset amb el qual podeu treballar és el “Heart Attack Analysis & Prediction dataset” (https://www.kaggle.com/datasets/rashikrahmanpritom/heart-attack-analysis-prediction-dataset).
Important: si escolliu un dataset diferent al proposat és important que aquest contingui una àmplia varietat de dades numèriques i categòriques per poder fer una anàlisi més rica i poder respondre a les diferents preguntes plantejades a l’enunciat de la pràctica. Seguint les principals etapes d’un projecte analític, les diferents tasques a realitzar (i justificar) són les següents:
Inicialment s’havia pensat aprofitar les dades recopilades durant la primera pràctica, relacionades amb dades de navegació de vaixells. No obstant, després de valorar-ho, no trobavem com poder utilitzar-lo de forma bastant gràfica i que ens permetés realitzar un estudi estadístic com el que es demana a l’enunciat, pel que finalment hem decidit utilitzar el conjunt de dades que es posa com a exemple a l’enunciat, ja que considerem que resulta molt interessant i pot donar més joc al estar relacionat amb temes de salut.
El conjunt de dades seleccionat “Heart Attack Analysis & Prediction dataset” proporciona informació sobre diferents factors que podrien estar relacionats amb malalties cardiovasculars, més concretament amb predir la probabilitat de patir un atac de cor.
Es tracta d’un conjunt de dades que ens permetrà realitzar un estudi en major profunditat de possibles problemes del cor, molt més centrat en els atacs de cor, al tenir una variable que està directament relacionada amb el resultat que es va obtenir durant l’estudi dels pacients i que ens permeti trobar un model per determinar la probabilitat de patir un atac de cor segons els valors dels diferents factors d’estudi.
Per poder facilitar aquesta anàlisi, el dataset conté el fitxer heart.csv amb les següents variables/atributs:
Per altra banda, hi ha disponible un segon fitxer, o2Saturation.csv, que conté múltiples observacions relacionades amb el nivell de saturació d’oxigen (una única variable). El fitxer conté moltíssimes més observacions (3586 en total) sense cap tipus d’idenfificador que ens permeti poder realitzar una integració de les dades d’ambdós fitxers per tenir un conjunt més complet.
Abans de començar, és necessari carregar el fitxer en un data frame que ens permeti treballar de forma còmoda amb les dades a través del codi en R:
# El joc de dades està conformat per un fitxer CSV, amb separació per comes i amb una capçalera, pel que carreguem el fitxer utilitzant la següent funció:
dfHeartAttack <- read.csv("../dataset/heart.csv", header=TRUE)
Una vegada carregades totes les dades, procedim a un primer anàlisi ràpid a partir del resum estadístic del data frame:
# Mostrem la informació bàsica del data frame per fer una comprovació a simple vista:
str(dfHeartAttack)
## 'data.frame': 303 obs. of 14 variables:
## $ age : int 63 37 41 56 57 57 56 44 52 57 ...
## $ sex : int 1 1 0 1 0 1 0 1 1 1 ...
## $ cp : int 3 2 1 1 0 0 1 1 2 2 ...
## $ trtbps : int 145 130 130 120 120 140 140 120 172 150 ...
## $ chol : int 233 250 204 236 354 192 294 263 199 168 ...
## $ fbs : int 1 0 0 0 0 0 0 0 1 0 ...
## $ restecg : int 0 1 0 1 1 1 0 1 1 1 ...
## $ thalachh: int 150 187 172 178 163 148 153 173 162 174 ...
## $ exng : int 0 0 0 0 1 0 0 0 0 0 ...
## $ oldpeak : num 2.3 3.5 1.4 0.8 0.6 0.4 1.3 0 0.5 1.6 ...
## $ slp : int 0 0 2 2 2 1 1 2 2 2 ...
## $ caa : int 0 0 0 0 0 0 0 0 0 0 ...
## $ thall : int 1 2 2 2 2 1 2 3 3 2 ...
## $ output : int 1 1 1 1 1 1 1 1 1 1 ...
# Mostrem també un resum estadístic de les variables del joc de dades:
skim(dfHeartAttack)
| Name | dfHeartAttack |
| Number of rows | 303 |
| Number of columns | 14 |
| _______________________ | |
| Column type frequency: | |
| numeric | 14 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| age | 0 | 1 | 54.37 | 9.08 | 29 | 47.5 | 55.0 | 61.0 | 77.0 | ▁▆▇▇▁ |
| sex | 0 | 1 | 0.68 | 0.47 | 0 | 0.0 | 1.0 | 1.0 | 1.0 | ▃▁▁▁▇ |
| cp | 0 | 1 | 0.97 | 1.03 | 0 | 0.0 | 1.0 | 2.0 | 3.0 | ▇▃▁▅▁ |
| trtbps | 0 | 1 | 131.62 | 17.54 | 94 | 120.0 | 130.0 | 140.0 | 200.0 | ▃▇▅▁▁ |
| chol | 0 | 1 | 246.26 | 51.83 | 126 | 211.0 | 240.0 | 274.5 | 564.0 | ▃▇▂▁▁ |
| fbs | 0 | 1 | 0.15 | 0.36 | 0 | 0.0 | 0.0 | 0.0 | 1.0 | ▇▁▁▁▂ |
| restecg | 0 | 1 | 0.53 | 0.53 | 0 | 0.0 | 1.0 | 1.0 | 2.0 | ▇▁▇▁▁ |
| thalachh | 0 | 1 | 149.65 | 22.91 | 71 | 133.5 | 153.0 | 166.0 | 202.0 | ▁▂▅▇▂ |
| exng | 0 | 1 | 0.33 | 0.47 | 0 | 0.0 | 0.0 | 1.0 | 1.0 | ▇▁▁▁▃ |
| oldpeak | 0 | 1 | 1.04 | 1.16 | 0 | 0.0 | 0.8 | 1.6 | 6.2 | ▇▂▁▁▁ |
| slp | 0 | 1 | 1.40 | 0.62 | 0 | 1.0 | 1.0 | 2.0 | 2.0 | ▁▁▇▁▇ |
| caa | 0 | 1 | 0.73 | 1.02 | 0 | 0.0 | 0.0 | 1.0 | 4.0 | ▇▃▂▁▁ |
| thall | 0 | 1 | 2.31 | 0.61 | 0 | 2.0 | 2.0 | 3.0 | 3.0 | ▁▁▁▇▆ |
| output | 0 | 1 | 0.54 | 0.50 | 0 | 0.0 | 1.0 | 1.0 | 1.0 | ▇▁▁▁▇ |
Veiem que tenim un joc de dades compost per un total de 303 observacions i 14 variables.
Realitzem un primer anàlisi exploratori que ens permeti entendre millor el conjunt de dades:
# Visualitzem la quantitat d'observacions que s'han identificat amb una probabilitat més alta i més baixa que es van identificar en referència als atacs de cor:
ggplot(dfHeartAttack,aes(factor(output, levels = c(0, 1), labels = c("Baixa", "Alta")), fill=factor(output, levels = c(0, 1), labels = c("Baixa", "Alta")))) +
geom_bar() +labs(x="Probabilitat", y="Pacients") +
guides(fill=guide_legend(title="")) +
scale_fill_manual(values=c("#008000","#800000")) +
ggtitle("Probabilitat d'atac de cor") +
theme(axis.text.x = element_text(hjust = 1))
# Calculem la proporció per poder tenir una idea de com es reparteixen les observacions:
table(dfHeartAttack$output)/length(dfHeartAttack$output)
##
## 0 1
## 0.4554455 0.5445545
Veiem que les dades estan molt repartides entre els pacients que tenen una probabilitat major i menor de patir un atac de cor segons els factors registrats.
Vegem la probabilitat de patir un atac de cor en base a les diferents variables:
plotbyAge <- ggplot(dfHeartAttack,aes(age, fill=factor(output, levels = c(0, 1), labels = c("Baixa", "Alta")))) +
geom_bar() +labs(x="Edat", y="Pacients") +
guides(fill=guide_legend(title="")) +
scale_fill_manual(values=c("#008000","#800000")) +
ggtitle("Probabilitat per edat") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plotbySex <- ggplot(dfHeartAttack,aes(factor(sex, levels = c(0, 1), labels = c("Femení", "Masculí")), fill=factor(output, levels = c(0, 1), labels = c("Baixa", "Alta")))) +
geom_bar() +labs(x="Gènere", y="Pacients") +
guides(fill=guide_legend(title="")) +
scale_fill_manual(values=c("#008000","#800000")) +
ggtitle("Probabilitat per gènere") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plotbyCP <- ggplot(dfHeartAttack,aes(factor(cp, levels = c(0, 1, 2, 3), labels = c("Angina típica", "Angina atípica", "No angina", "Assimptomàtic")), fill=factor(output, levels = c(0, 1), labels = c("Baixa", "Alta")))) +
geom_bar() +labs(x="Tipus de dolor", y="Pacients") +
guides(fill=guide_legend(title="")) +
scale_fill_manual(values=c("#008000","#800000")) +
ggtitle("Probabilitat per dolor toràcic") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plotbyExng <- ggplot(dfHeartAttack,aes(factor(exng, levels = c(0, 1), labels = c("No", "Sí")), fill=factor(output, levels = c(0, 1), labels = c("Baixa", "Alta")))) +
geom_bar() +labs(x="Angina provocada per exercici", y="Pacients") +
guides(fill=guide_legend(title="")) +
scale_fill_manual(values=c("#008000","#800000")) +
ggtitle("Probabilitat per angina provocada per exercici") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plotbyFBS <- ggplot(dfHeartAttack,aes(factor(fbs, levels = c(0, 1), labels = c("Normal", "Alt")), fill=factor(output, levels = c(0, 1), labels = c("Baixa", "Alta")))) +
geom_bar() +labs(x="Nivell de sucre", y="Pacients") +
guides(fill=guide_legend(title="")) +
scale_fill_manual(values=c("#008000","#800000")) +
ggtitle("Probabilitat per nivell de sucre") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plotbyTRTBPS <- ggplot(dfHeartAttack,aes(trtbps, fill=factor(output, levels = c(0, 1), labels = c("Baixa", "Alta")))) +
geom_bar() +labs(x="Pressió arterial", y="Pacients") +
guides(fill=guide_legend(title="")) +
scale_fill_manual(values=c("#008000","#800000")) +
ggtitle("Probabilitat per pressió arterial en repòs") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plotbyChol <- ggplot(dfHeartAttack,aes(chol, fill=factor(output, levels = c(0, 1), labels = c("Baixa", "Alta")))) +
geom_bar() +labs(x="Colesterol", y="Pacients") +
guides(fill=guide_legend(title="")) +
scale_fill_manual(values=c("#008000","#800000")) +
ggtitle("Probabilitat per colesterol") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plotbyRestecg <- ggplot(dfHeartAttack,aes(factor(restecg, levels = c(0, 1, 2), labels = c("Normal", "Anomalies", "Hipertròfia")), fill=factor(output, levels = c(0, 1), labels = c("Baixa", "Alta")))) +
geom_bar() +labs(x="Electrocardiograma", y="Pacients") +
guides(fill=guide_legend(title="")) +
scale_fill_manual(values=c("#008000","#800000")) +
ggtitle("Probabilitat per resultat de l'electrocardiograma") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plotbyOldpeak <- ggplot(dfHeartAttack,aes(oldpeak, fill=factor(output, levels = c(0, 1), labels = c("Baixa", "Alta")))) +
geom_bar() +labs(x="Canvis electrocardiograma", y="Pacients") +
guides(fill=guide_legend(title="")) +
scale_fill_manual(values=c("#008000","#800000")) +
ggtitle("Probabilitat per canvis en l'electrocardiograma") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plotbySLP <- ggplot(dfHeartAttack,aes(factor(slp, levels = c(0, 1, 2), labels = c("Plana", "Ascendent", "Descendent")), fill=factor(output, levels = c(0, 1), labels = c("Baixa", "Alta")))) +
geom_bar() +labs(x="Pendent", y="Pacients") +
guides(fill=guide_legend(title="")) +
scale_fill_manual(values=c("#008000","#800000")) +
ggtitle("Probabilitat per patró de canvi") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plotbyThalachh <- ggplot(dfHeartAttack,aes(thalachh, fill=factor(output, levels = c(0, 1), labels = c("Baixa", "Alta")))) +
geom_bar() +labs(x="Freqüència", y="Pacients") +
guides(fill=guide_legend(title="")) +
scale_fill_manual(values=c("#008000","#800000")) +
ggtitle("Probabilitat per freqüència màxima registrada") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plotbyCAA <- ggplot(dfHeartAttack,aes(caa, fill=factor(output, levels = c(0, 1), labels = c("Baixa", "Alta")))) +
geom_bar() +labs(x="Vasos", y="Pacients") +
guides(fill=guide_legend(title="")) +
scale_fill_manual(values=c("#008000","#800000")) +
ggtitle("Probabilitat per vasos sanguinis obstruïts") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plotbyThall <- ggplot(dfHeartAttack,aes(factor(thall, levels = c(0, 1, 2, 3), labels = c("Unknown", "Sense", "Defecte fix", "Defecte reversible")), fill=factor(output, levels = c(0, 1), labels = c("Baixa", "Alta")))) +
geom_bar() +labs(x="Indicis", y="Pacients") +
guides(fill=guide_legend(title="")) +
scale_fill_manual(values=c("#008000","#800000")) +
ggtitle("Probabilitat per talassèmia") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
grid.newpage()
grid.arrange(plotbyAge, plotbySex, plotbyCP, plotbyExng, plotbyFBS, plotbyTRTBPS, plotbyChol, plotbyRestecg, plotbyOldpeak, plotbySLP, plotbyThalachh, plotbyCAA, plotbyThall, ncol=2)
Finalment, creem un gràfic on es mostrin les possibles relacions existentes entre les variables per veure, de forma ràpida, si podem veure algunes correlacions de forma directa:
# Generem un gràfic que compara parells de variables entre elles, el que ens permet
ggpairs(data = dfHeartAttack, columns = 1:13)
Observant els diversos factors que es troben dins el conjunt de dades i tots els gràfics creats, sembla que la variància entre els diferents valors, així com la, a priori, poca correlació entre les variables i la seva relació amb la probabilitat de patir un atac de cor podria ser interessant mantenir tots els atributs que tenim per poder utilitzar-los durant la resta de l’anàlisi.
La integració consisteix en la combinació de dades de diferents fonts, per tal de crear una estructura de dades coherent. En el cas d’estudi, aquesta integració ja està realitzada i tenim, de cada observació, totes les variables en columnes.
La selecció consisteix en filtrar o seleccionar les dades d’interès. A partir de l’anàlisi exploratori que hem realitzat consideram vàlides totes les dades i no cal filtrar ni reduir la dimensionalitat del dataset.
Comprobem si existeixen valors nul·ls en les dades, tot i que, aparentment, al revisar els resultats del resum estadístic anterior sembla que no existeixi cap valor d’aquest tipus:
# Revisem el total de valors NA i blancs que hi ha a cada variable:
colSums(is.na(dfHeartAttack))
## age sex cp trtbps chol fbs restecg thalachh
## 0 0 0 0 0 0 0 0
## exng oldpeak slp caa thall output
## 0 0 0 0 0 0
colSums(dfHeartAttack == "")
## age sex cp trtbps chol fbs restecg thalachh
## 0 0 0 0 0 0 0 0
## exng oldpeak slp caa thall output
## 0 0 0 0 0 0
Es confirma que no existeixen valors d’aquest tipus en el conjunt de dades, pel que no serà necessari realitzar cap acció al respecte. En cas de tenir valors d’aquest tipus hauriem de pensar si treure’ls o bé realitzar una aproximació del possible valor a partir de la resta de valors de la variable en qüestió.
Respecte a les variables categòriques, podem comprovar en el resum estadístic que no hi ha cap valor fora del rang vàlid; per tant no tenen valors atípics.
Respecte a les variables numèriques, utilitzarem els diagrames de caixa per poder veure ràpidament si existeix algun valor atípic o outlier en el joc de dades:
# Calculem els diferents gràfics de caixa:
bpAge <- ggplot(data = dfHeartAttack) +
geom_boxplot(aes(y = age, fill = "Edat")) +
scale_fill_manual(values = "#008080") +
ggtitle("Boxplot per edat")
bpTRTBPS <- ggplot(data = dfHeartAttack) +
geom_boxplot(aes(y = trtbps, fill = "Pres.")) +
scale_fill_manual(values = "#00A595") +
ggtitle("Boxplot per pressió arterial en repòs")
bpChol <- ggplot(data = dfHeartAttack) +
geom_boxplot(aes(y = chol, fill = "Col.")) +
scale_fill_manual(values = "#00B080") +
ggtitle("Boxplot per nivell de colesterol")
bpOldpeak <- ggplot(data = dfHeartAttack) +
geom_boxplot(aes(y = oldpeak, fill = "Canv.")) +
scale_fill_manual(values = "#00C0B0") +
ggtitle("Boxplot per canvis en l'electrocardiograma")
bpThalachh <- ggplot(data = dfHeartAttack) +
geom_boxplot(aes(y = thalachh, fill = "Freq.")) +
scale_fill_manual(values = "#00D0A0") +
ggtitle("Boxplot per freqüència màxima registrada")
# Mostrem els diferents gràfics de forma agrupada:
grid.newpage()
grid.arrange(bpAge, bpTRTBPS, bpChol, bpOldpeak, bpThalachh, ncol=3)
S’observen valors extrems a partir dels gràfics per algunes de les variables. No obstant, veiem que són valors que estan dintre dels paràmetres que es poden considerar com a vàlids:
En resum, no s’han trobat valors anòmals que puguin considerar-se fora dels valors possibles, tot i que sí hi ha alguns valors extrems degut a les condicions físiques i/o de salut dels diferents pacients. Possiblement aquests valors poden ser importants a l’hora de fer estimacions. No procedeix eliminar cap valor atípic.
La normalització de les dades ens permet obtenir valors en escales que permetin comparar la magnitud de forma similar entre els diferents rangs de valors que tenen les variables.
Per altra banda, la discretització ens permet agrupar observacions numèriques per tenir noves categories que puguin resultar útils durant l’anàlisi (per exemple en algorismes de classificació d’arbre):
# Seleccionem la llavor per l'algoritme randomitzador:
set.seed(1234)
# data frame Heart Attack Normalitzat (dfHAN)
dfHAN <- dfHeartAttack %>%
mutate( # Discretització
age.d = discretize(age, "cluster", breaks = 5),
trtbps.d = discretize(trtbps, "cluster", breaks = 3),
chol.d = discretize(chol, "cluster", breaks = 5),
thalachh.d = discretize(thalachh, "cluster", breaks = 3),
# Normalització
age = scale( age, center=TRUE, scale=TRUE )[,1],
sex = factor( ifelse( sex==1, "Masculí",
ifelse( sex==0, "Femení", "?"))),
cp = factor( ifelse( cp==0, "Angina típica",
ifelse( cp==1, "Angina atípica",
ifelse( cp==2, "Dolor no relacionat amb angina",
ifelse( cp==3, "Sense dolor toràcic", "?"))))),
trtbps = scale( trtbps, center=TRUE, scale=TRUE )[,1],
chol = scale( chol, center=TRUE, scale=TRUE )[,1],
fbs = factor( ifelse( fbs==0, "Normal",
ifelse( fbs==1, "Alt", "?"))),
restecg = factor( ifelse( restecg==0, "Normal",
ifelse( restecg==1, "Anomalies ST-T",
ifelse( restecg==2, "Hipertrofia ventricular", "?")))),
thalachh = scale( thalachh, center=TRUE, scale=TRUE )[,1],
exng = factor( ifelse( exng==0, "No exercici",
ifelse( exng==1, "Si exercici", "?"))),
oldpeak = scale( oldpeak, center=TRUE, scale=TRUE )[,1],
slp = factor( ifelse( slp==0, "Pendent plana",
ifelse( slp==1, "Pendent ascendent",
ifelse( slp==2, "Pendent descendent", "?")))),
caa = factor( ifelse( caa==0, "Sense obstrucció",
ifelse( caa==1, "Obstrucció en un vas",
ifelse( caa==2, "Obstrucció de dos vasos",
ifelse( caa==3, "Obstrucció de tres vasos",
ifelse( caa==4, "Obstrucció de quatre vasos", "?")))))),
thall = factor( ifelse( thall==1, "No hi ha antecedents",
ifelse( thall==2, "Presència defecte fix",
ifelse( thall==3, "Presència de defecte reversible", "?")))),
resultat = factor( ifelse( output==0, "Atac probable",
ifelse( output==1, "Atac poc probable", "?"))))
dfHANnum <- dfHAN %>%
select(age, trtbps, chol, thalachh, oldpeak )
dfHANcat <- dfHAN %>%
select(age.d, trtbps.d, chol.d, thalachh.d, sex, cp, fbs, restecg, exng, slp, caa, thall, resultat)
skim(dfHAN)
| Name | dfHAN |
| Number of rows | 303 |
| Number of columns | 19 |
| _______________________ | |
| Column type frequency: | |
| factor | 13 |
| numeric | 6 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| sex | 0 | 1 | FALSE | 2 | Mas: 207, Fem: 96 |
| cp | 0 | 1 | FALSE | 4 | Ang: 143, Dol: 87, Ang: 50, Sen: 23 |
| fbs | 0 | 1 | FALSE | 2 | Nor: 258, Alt: 45 |
| restecg | 0 | 1 | FALSE | 3 | Ano: 152, Nor: 147, Hip: 4 |
| exng | 0 | 1 | FALSE | 2 | No : 204, Si : 99 |
| slp | 0 | 1 | FALSE | 3 | Pen: 142, Pen: 140, Pen: 21 |
| caa | 0 | 1 | FALSE | 5 | Sen: 175, Obs: 65, Obs: 38, Obs: 20 |
| thall | 0 | 1 | FALSE | 4 | Pre: 166, Pre: 117, No : 18, ?: 2 |
| age.d | 0 | 1 | FALSE | 5 | [53: 104, [61: 71, [46: 57, [40: 52 |
| trtbps.d | 0 | 1 | FALSE | 3 | [11: 153, [13: 97, [94: 53 |
| chol.d | 0 | 1 | FALSE | 5 | [23: 101, [19: 100, [28: 65, [12: 32 |
| thalachh.d | 0 | 1 | FALSE | 3 | [13: 141, [16: 85, [71: 77 |
| resultat | 0 | 1 | FALSE | 2 | Ata: 165, Ata: 138 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| age | 0 | 1 | 0.00 | 1.0 | -2.79 | -0.76 | 0.07 | 0.73 | 2.49 | ▁▆▇▇▁ |
| trtbps | 0 | 1 | 0.00 | 1.0 | -2.15 | -0.66 | -0.09 | 0.48 | 3.90 | ▃▇▅▁▁ |
| chol | 0 | 1 | 0.00 | 1.0 | -2.32 | -0.68 | -0.12 | 0.54 | 6.13 | ▃▇▂▁▁ |
| thalachh | 0 | 1 | 0.00 | 1.0 | -3.43 | -0.70 | 0.15 | 0.71 | 2.29 | ▁▂▅▇▂ |
| oldpeak | 0 | 1 | 0.00 | 1.0 | -0.90 | -0.90 | -0.21 | 0.48 | 4.44 | ▇▂▁▁▁ |
| output | 0 | 1 | 0.54 | 0.5 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00 | ▇▁▁▁▇ |
Una possible forma de reduir-ne la dimensionalitat és considerar els components principals de les variables numériques:
pca.acc <- prcomp(dfHANnum, scale. = TRUE)
summary( pca.acc )
## Importance of components:
## PC1 PC2 PC3 PC4 PC5
## Standard deviation 1.3441 1.0380 0.9399 0.8713 0.68799
## Proportion of Variance 0.3613 0.2155 0.1767 0.1518 0.09467
## Cumulative Proportion 0.3613 0.5768 0.7535 0.9053 1.00000
Podem veure que partim de cinc variables i necessitam quatre per arribar a descriure el 90% de la variabilitat total, per la qual cosa no suposa una gran reducció de dimensionalitat.
Per facilitar els càlculs futurs, es procedirà a un anàlisi de la normalitat dels valors numèrics del conjunt de dades, el que ens permetrà saber si és possible aplicar certs tests més endavant.
Comencem amb un test de Shapiro-Wilk de normalitat:
# Realitzem el test de normalitat Shapiro-Wilk a les diferents variables numériques
apply( dfHANnum, 2, shapiro.test )
## $age
##
## Shapiro-Wilk normality test
##
## data: newX[, i]
## W = 0.98637, p-value = 0.005798
##
##
## $trtbps
##
## Shapiro-Wilk normality test
##
## data: newX[, i]
## W = 0.96592, p-value = 1.458e-06
##
##
## $chol
##
## Shapiro-Wilk normality test
##
## data: newX[, i]
## W = 0.94688, p-value = 5.365e-09
##
##
## $thalachh
##
## Shapiro-Wilk normality test
##
## data: newX[, i]
## W = 0.97632, p-value = 6.621e-05
##
##
## $oldpeak
##
## Shapiro-Wilk normality test
##
## data: newX[, i]
## W = 0.84418, p-value < 2.2e-16
Tots els test mostren que les variables no segueixen una distribució totalment normal. No obstant, pel Teorema del Límit Central, al tractar-se d’una quantitat elevada d’observacions, es pot considerar que tendeixen a una distribució normal. Per corroborar aquest punt, generem els gràfics QQ per veure si la tendència és normal o no:
# Creem els diferents gràfics QQ per comprovar si les dades de les variables segueixen o s'aproximen a una distribució normal:
qqAge <- ggplot(data = dfHAN, aes(sample = age)) +
geom_qq() +
geom_qq_line() +
labs(x = "Valors teòrics", y = "Valors observats") +
ggtitle("Gràfic QQnorm: 'age'")
qqTRTBPS <- ggplot(data = dfHAN, aes(sample = trtbps)) +
geom_qq() +
geom_qq_line() +
labs(x = "Valors teòrics", y = "Valors observats") +
ggtitle("Gràfic QQnorm: 'trtbps'")
qqChol <- ggplot(data = dfHAN, aes(sample = chol)) +
geom_qq() +
geom_qq_line() +
labs(x = "Valors teòrics", y = "Valors observats") +
ggtitle("Gràfic QQnorm: 'chol'")
qqOldpeak <- ggplot(data = dfHAN, aes(sample = oldpeak)) +
geom_qq() +
geom_qq_line() +
labs(x = "Valors teòrics", y = "Valors observats") +
ggtitle("Gràfic QQnorm: 'oldpeak'")
qqThalachh <- ggplot(data = dfHAN, aes(sample = thalachh)) +
geom_qq() +
geom_qq_line() +
labs(x = "Valors teòrics", y = "Valors observats") +
ggtitle("Gràfic QQnorm: 'thalachh'")
# Mostrem els diferents gràfics de forma agrupada:
grid.newpage()
grid.arrange(qqAge, qqTRTBPS, qqChol, qqOldpeak, qqThalachh, ncol=3)
Amb els resultats obtinguts, podem observar que les variables numèriques, tot i no tenir un a distribució normal, sí semblen tendir a aquesta. L’únic cas que podriem considerar que no és tant així seria per la variable oldpeak, ja que és la única que no té una distribució tant propera a la línia que podriem considerar que defineix la normalitat.
Tot i que ja s’ha vist anteriorment que semblava que no hi havia correlacions molt fortes entre les diferents variables, procedim a analitzar-ho més a fons:
# Generam un gràfic de correlació que ens permeti analitzar-ho fàcilment (només aquelles variables que no siguin factors):
res<-cor(dfHeartAttack[, 1:14])
corrplot(res, method="color", tl.col="black", tl.srt=30, order = "AOE", number.cex=0.75, sig.level = 0.01, addCoef.col = "black")
La més forta d’aquestes correspon amb una correlació negativa entre el patró de canvi de l’electrocardiograma en una situació d’estrés (variable slp) i el canvi en el segment ST de l’electrocardiograma (variable oldpeak). Mentre que la més forta amb valor positiu correspon a ula correlació entre el tipus de dolor toràcic (variable cp) i la probabilitat d’un atac de cor (variable output).
Cal tenir en compte que aquestes correlacions no són molt fortes i el signe ens indica que quan una augmenta o disminueix, la variable correlacionada actúa de la mateixa manera en certa proporció segons la potència d’aquesta correlació.
Podem fer la correlació únicament sobre les variables numériques:
# Generam un gràfic de correlació que ens permeti analitzar-ho fàcilment (només aquelles variables que no siguin factors):
res <- cor( dfHANnum )
corrplot(res, method="color", tl.col="black", tl.srt=30, order = "AOE", number.cex=0.75, sig.level = 0.01, addCoef.col = "black")
Amb els resultats obtinguts, no hi ha cap evidència que, per les correlacions existents, es pugui disminuïr el nombre de variables.
La regresió logística és un model que intentar predir el resultat
(output) a partir de la resta de variabels mitjançant un
model de regressió estadística.
regressio <- lm(output ~ age + sex + cp + trtbps + chol + fbs + restecg + thalachh + exng + oldpeak + slp + caa + thall, data = dfHAN)
summary(regressio)
##
## Call:
## lm(formula = output ~ age + sex + cp + trtbps + chol + fbs +
## restecg + thalachh + exng + oldpeak + slp + caa + thall,
## data = dfHAN)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.01468 -0.18519 0.03069 0.22575 1.02624
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.30215 0.25796 1.171 0.242466
## age 0.02464 0.02406 1.024 0.306638
## sexMasculí -0.16407 0.04823 -3.402 0.000766 ***
## cpAngina típica -0.16886 0.06397 -2.640 0.008767 **
## cpDolor no relacionat amb angina 0.05637 0.06102 0.924 0.356389
## cpSense dolor toràcic 0.09685 0.08875 1.091 0.276084
## trtbps -0.04104 0.02137 -1.920 0.055844 .
## chol -0.01653 0.02084 -0.793 0.428416
## fbsNormal -0.03335 0.05749 -0.580 0.562362
## restecgHipertrofia ventricular -0.13405 0.17685 -0.758 0.449098
## restecgNormal -0.04699 0.04083 -1.151 0.250823
## thalachh 0.04285 0.02563 1.672 0.095634 .
## exngSi exercici -0.09352 0.04999 -1.871 0.062437 .
## oldpeak -0.04828 0.02674 -1.806 0.072023 .
## slpPendent descendent 0.13900 0.04953 2.806 0.005364 **
## slpPendent plana 0.06540 0.08593 0.761 0.447281
## caaObstrucció de quatre vasos 0.38693 0.16945 2.284 0.023149 *
## caaObstrucció de tres vasos 0.04660 0.09558 0.488 0.626212
## caaObstrucció en un vas 0.07225 0.07324 0.987 0.324711
## caaSense obstrucció 0.34261 0.06825 5.020 9.2e-07 ***
## thallNo hi ha antecedents 0.22837 0.25437 0.898 0.370068
## thallPresència de defecte reversible 0.07849 0.24369 0.322 0.747623
## thallPresència defecte fix 0.28800 0.24183 1.191 0.234694
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.335 on 280 degrees of freedom
## Multiple R-squared: 0.5818, Adjusted R-squared: 0.549
## F-statistic: 17.71 on 22 and 280 DF, p-value: < 2.2e-16
Podem veure que amb la regressió múltiple \(R^2=0.5489541\), que ens indica que, amb
aquest model, tot el conjunt de variables disponibles expliquen el 54.9%
de la variable objetivo output.
Ens demanam si les dones tenen major probabilitat de tenir un ataca de cor que els dones. Podem veure la gràfica comparativa per gènere:
plotbySex
En la mostra hi ha un percentatge major de dones amb alta probabilitat de tenir un atac de cor que d’homes; però ens demanam si aquesta diferència és significativa ja que podria ser fruït de l’atzar de la mostra.
La hipòtesi nul·la, o \(H_{0}\) és: La probabilitat de tenir un atac de cor de les dones és igual a la probabilitat de tenir un atac de cor dels homes
La hipòtesi alternativa, o \(H_{1}\) és: La probabilitat de tenir un atac de cor de les dones és superior a la probabilitat de tenir un atac de cor dels homes
Es tracta d’un test de contraste unilateral sobre dues mostres (les probabilitats de tenir atac de cor de les dones i les dels homes) en relació a les seves probabilitats de tenir un atac de cor.
Per determinar el test a aplicar, és pertinent aplicar el teorema del límit central que estableix que el contrast d’hipòtesi sobre una mitjana d’una mostra s’aproxima a una distribució normal encara que la població original no segueixi una distribució normal, sempre que la mida de la mostra sigui suficientment gran (superior a 30). Això ja ho hem vist en l’apartat del Test de normalitat.
Perquè es puguin donar les condicions para aplicar el teorema del límit central, el tamany de les dues mostres ha de ser superior a 30; vegem-ho:
mostra_dones <- dfHAN$output[dfHAN$sex=="Femení"]
mostra_homes <- dfHAN$output[dfHAN$sex=="Masculí"]
print( paste("La mostra de les dones és de", length(mostra_dones), "observacions") )
## [1] "La mostra de les dones és de 96 observacions"
print( paste("La mostra dels homes és de", length(mostra_homes), "observacions") )
## [1] "La mostra dels homes és de 207 observacions"
Podem veure que les dues mostres compleixen les hipòtesis del teorema del límit central i, per tant podem assumir que segueixen lleis normals.
Les variances de les poblacions (la població de les dones i dels homes) no les coneixem. Només tenim, òbviament, les respectives variances mostrals. A continuació, necessitam saber si les variances de les dues poblacions són iguales o no, Per això cal fer un test d’homoscedasticitat:
Assumim (per l’argument anterior) que les dues mostres correponen a dues poblacions normals independents \(N(\mu_{1},\sigma_{1})\) i \(N(\mu_{2},\sigma_{2})\). Aleshores la variable aleatòria següent segueix una distribució F d’Snedecor amb \(n_1-1\) i \(n_2-1\) graus de llibertat: \[ F=\frac{\frac{s_{1}^{2}}{\sigma_{1}^{2}}}{\frac{s_{2}^{2}}{\sigma_{2}^{2}}} \]
On \(s_1\) i \(s_2\) són les desviacions estàndards mostrals i \(\sigma_{1}\) i \(\sigma_{2}\) les desviacions poblacionals.
Sota la hipòtesi nul·la \(H_{0}: \sigma_{1}^{2}=\sigma_{2}^{2}\), el test estadístic és: \[ f_{obs}=\frac{s_{1}^{2}}{s_{2}^{2}}\sim F_{n_{1}-1,n_{2}-1} \] On F és una distribució d’Snedecor amb \(n_1-1\) i \(n_2-1\) graus de llibertat.
La hipòtesi alternativa \(H_{1}: \sigma_{1}^{2}\neq \sigma_{2}^{2}\)
Facem el càlcul en R:
var.test( mostra_dones, mostra_homes)
##
## F test to compare two variances
##
## data: mostra_dones and mostra_homes
## F = 0.76208, num df = 95, denom df = 206, p-value = 0.1343
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.5455293 1.0885394
## sample estimates:
## ratio of variances
## 0.7620767
Podem veure que el valor de l’estadístic cau en la zona d’acceptació de l’hipòtesi nul·la. De fet, obtenim un valor de p superior a \(\alpha\) i, per tant, no podem rebutjar la hipòtesi nul·la. Per tant, les variances poblacionals són essencialment iguales.
En resum, el test a aplicar és el corresponent a dues mostres independents corresponents a poblacions que se poden aproximar a una llei normal, sobre la mitjana amb variances desconegudes però iguales.
El test estadístic de la diferència de les dues mitjanes segueix una distribució t d’Student amb \(n_{1}+n_{2}-2\) graus de llibertat: \[ t=\frac{\bar{x}_{1}-\bar{x}_{2}}{S\sqrt{\frac{1}{n_{1}}+\frac{1}{n_{2}}}}\sim t_{n_{1}+n_{2}-2} \]
El valor S és la desviació típica comuna que es calcula com: \[ S=\sqrt{\frac{(n_1-1)S_1^2+(n_2-1)S_2^2}{n_1+n_2-2}} \] On \(s_1^2\) i \(s_2^2\) són les variances mostrals estimades.
Sota la hipòtesi nul·la \(H_0\): \(P(t_{\alpha/2}\le t_{obs} \le t_{1-\alpha/2})=1-\alpha\)
Per tant, la zona d’acceptació és \([t_{\alpha/2},t_{1-\alpha/2}]\).
Facem el càlcul en R:
alfa <- 0.05
# Tamany de les mostres
nd <- length(mostra_dones)
nh <- length(mostra_homes)
# Mitjanes de les mostres
md <- mean(mostra_dones)
mh <- mean(mostra_homes)
# Desviació estàndard de les mostres
sd <- sd(mostra_dones)
sh <- sd(mostra_homes)
# Valor de l'estadístic
S <- sqrt( ( (nd-1)*sd^2 + (nh-1)*sh^2 ) / (nd+nh-2) )
tobs <- (md-mh)/(S*sqrt( 1/nd + 1/nh ))
# Valors crítics
tcrit.L <- qt( alfa/2, df=nd+nh-2)
tcrit.U <- qt( alfa/2, df=nd+nh-2, lower.tail=FALSE)
# Valor de p
valor_p <- pt( abs(tobs), df=nd+nh-2, lower.tail=FALSE)*2
print(paste("Valor observat:", round(tobs,2)))
## [1] "Valor observat: 5.08"
print(paste("Interval d'acceptació: [",round(tcrit.L, 2),",", round(tcrit.U, 2), "]"))
## [1] "Interval d'acceptació: [ -1.97 , 1.97 ]"
print(paste("Valor de p:", round(valor_p, 2)))
## [1] "Valor de p: 0"
Podem veure que l’observació de l’estadístic considerat és 5.08, valor que cau fora la zona d’acceptació: [-1.97, 1.97]. Per tant podem rebutjar la hipòtesi nul·la amb un nivell de confiança del 95%.
El valor de p (probabilitat de l’error que s’estaria cometent si es rebutja la hipòtesi nul·la essent aquesta certa) és de 0. Un valor inferior a \(\alpha\).
Com a conseqüència, hem de rebutjar la hipòtesi nul·la a favor de l’alternativa: La probabilitat de tenir un atac de cor de les dones és superior a la probabilitat de tenir un atac de cor dels homes
Anem a generar un model no supervisat basat en l’algorisme k_means. El paràmetre fundamental és el valor de k, que coincideix amb el nombre de grups que volem trobar. En realitat en el nostre problema volem classificar les observacions en dos grups: els que tenen risc d’atac de cor i els que no; per tant, el valor que ens interessa és k=2; però abans analitzarem quin és el valor de k amb el que podem obtenir uns clústers més diferenciats:
Estimació del millor valor de k pel mètode Calinski-Harabasz
fit_ch <- kmeansruns(dfHANnum, krange = 1:10, criterion = "ch")
fit_asw <- kmeansruns(dfHANnum, krange = 1:10, criterion = "asw")
fit_ch$bestk
## [1] 2
fit_asw$bestk
## [1] 2
plot(1:10,fit_ch$crit,type="o",col="blue",pch=0,xlab="Número de clústers",ylab="Criteri Calinski-Harabasz")
plot(1:10,fit_asw$crit,type="o",col="blue",pch=0,xlab="Número de clústers",ylab="Criteri silueta mitja")
Podem veure que el millor valor de k és 2, per tant, facem una clusterització en 2 grups i analitzem la qualitat. Per poder visualitzar el resultat, ho projectarem sobre les dues primeres components principals:
d <- daisy(dfHANnum)
millor_qualitat <- 0
for (j in 1:50) {
fit <- kmeans(dfHANnum, 2)
y_cluster <- fit$cluster
sk <- silhouette(y_cluster, d)
qualitat <- mean(sk[,3])
if (qualitat>millor_qualitat) {
millor_qualitat <- qualitat
millor_fit <- fit
}
}
fit2 <- millor_fit
y_cluster2 <- fit2$cluster
clusplot(dfHANnum, fit2$cluster,
color=TRUE, shade=TRUE, labels=5, lines=0)
sk2 <- silhouette(y_cluster2, d)
qualitat2 <- mean(sk2[,3])
qualitat2
## [1] 0.2435173
La qualitat de la clusterització segons la silueta és 0.2435173; un valor molt baix. Visualment no s’aprecia una clara separació dels dos grups.
Comparem els dos grups trobats amb la variable
output:
resultat1 <- matrix( rep(0,4), nrow=2, ncol = 2,
dimnames=list( c("Grup 1", "Grup 2"),
c("Atac probable", "Atac poc probable")))
for (i in 1:nrow(dfHANnum)) {
resultat1[fit2$cluster[i],dfHAN$output[i]+1] <- resultat1[fit2$cluster[i],dfHAN$output[i]+1] + 1
}
total <- sum(resultat1)
opc1 <- round(100 * (resultat1["Grup 1", "Atac probable"] + resultat1["Grup 2", "Atac poc probable"] ) / total, 1)
opc2 <- round( 100 * (resultat1["Grup 1", "Atac poc probable"] + resultat1["Grup 2", "Atac probable"] ) / total, 1)
millor_opc <- max( opc1, opc2)
resultat1
## Atac probable Atac poc probable
## Grup 1 50 128
## Grup 2 88 37
print(paste("Si Grup 1 = Atac probable llavors, percentatge encert:", opc1 ))
## [1] "Si Grup 1 = Atac probable llavors, percentatge encert: 28.7"
print(paste("Si Grup 1 = Atac poc probable llavors, percentatge encert:", opc2))
## [1] "Si Grup 1 = Atac poc probable llavors, percentatge encert: 71.3"
Podem veure que la millor opció d’interpretació dels resultats dels dos grups ens dóna un percentatge d’encert de 71.3%.
Conclusió:
Els dos grups trobats considerats com a clústers no estan gens separats en distància. De fet si calculam la mesura de la silueta sobre la classificació coneguda, obtenim una qualitat molt baixa de 0.2393498. L’algorismes de clusterització només encerta en un 71.3% dels casos amb la classificació real.
Ara aplicarem un model supervisat basat en la generació de regles de classificació a partir d’un arbre de decisions. Farem servir les dades categòriques amb la discretització de les variables numériques; però abans hem de separar el data set en dos: un per entrenar l’algorisme i l’altre per validar-ho. Farem servir la proporció habitual de 2/3 per entrenament i 1/3 per test.
indexes = 1:(nrow(dfHAN)/3) * 3 - 2
train <- dfHANcat[indexes,]
test <- dfHANcat[-indexes,]
testY <- dfHAN[-indexes,14]
Una vegada feta la separació aleatòria de les mostres, convé realitzar una mínima anàlisi de dades per a assegurar-nos de no obtenir classificadors esbiaixats pels valors que conté cada mostra. En aquest cas, verificarem que la proporció del prèstecs bons i dolents és més o menys constant en els dos conjunts:
bons_tr <- ( train %>%
filter( resultat == "Atac probable") %>%
nrow() )
num_tr <- train %>% nrow()
percen_tr <- round((bons_tr / num_tr)*100, digits=1)
bons_ts <- ( test %>%
filter( resultat == "Atac probable") %>%
nrow() )
num_ts <- test %>% nrow()
percen_ts <- round((bons_ts / num_ts)*100, digits=1)
Es pot veure que els percentatges de probabilitat d’atacs de cor del conjunt d’entrenament i del conjunt de test són semblants (45.5% i 45.5% respectivament). Per tant, consideram correctes els jocs d’entrenament i el de test.
Anem a generar el model de regles de decisió per determinar si la probabilitat d’atac de cor és o no alta. Ho farem amb l¡’algorisme random Forest, tècnica en què genera diversos classificadors amb els seus arbres de decisió per, finalment, registrar-ne tots els resultats i determinar-ne la classe final.
rf <- randomForest(resultat ~ ., data=train)
print(rf)
##
## Call:
## randomForest(formula = resultat ~ ., data = train)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 3
##
## OOB estimate of error rate: 13.86%
## Confusion matrix:
## Atac poc probable Atac probable class.error
## Atac poc probable 50 5 0.09090909
## Atac probable 9 37 0.19565217
Validem el model amb el joc de test:
pred = predict(rf, newdata=test[1:12])
cm = table(test[,13], pred)
precisio_model <- round( 100 * sum(diag(cm)) / sum(cm), digits = 2 )
print(sprintf("La precisión del árbol es: %.1f %%", precisio_model))
## [1] "La precisión del árbol es: 80.2 %"
Podem veure que obtenim un model de decisió d’arbre amb una taxa de predicció del (80.2%).
Vegem la taula de confusió:
CrossTable(test[,13], pred, prop.chisq = FALSE, prop.c = FALSE, prop.r =TRUE,dnn = c('Reality', 'Prediction'))
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Row Total |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 202
##
##
## | Prediction
## Reality | Atac poc probable | Atac probable | Row Total |
## ------------------|-------------------|-------------------|-------------------|
## Atac poc probable | 92 | 18 | 110 |
## | 0.836 | 0.164 | 0.545 |
## | 0.455 | 0.089 | |
## ------------------|-------------------|-------------------|-------------------|
## Atac probable | 22 | 70 | 92 |
## | 0.239 | 0.761 | 0.455 |
## | 0.109 | 0.347 | |
## ------------------|-------------------|-------------------|-------------------|
## Column Total | 114 | 88 | 202 |
## ------------------|-------------------|-------------------|-------------------|
##
##
A partir del fitxer de dades, hem fet unes tasques de preprocessament:
Amb les dades preprocessades, hem fet una sèrie d’anàlisi:
En definitiva, amb aquesta pràctica hem seguit les fases del cicle de vida de les dades: captura, emmagazematge, prepocessat, anàlisi, visualització i publicació.